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Observations of galaxies over large distances reveal the possibility of a fractal distribution of their 
positions. The source of fractal behavior is the lack of a length scale in the two body gravitational 
interaction. However, even with new, larger, sample sizes from recent surveys, it is difficult to extract 
information concerning fractal properties with confidence. Similarly, simulations with a billion par- 
ticles only provide a thousand particles per dimension, far too small for accurate conclusions. With 
one dimensional models these limitations can be overcome by carrying out simulations with on the 
order of a quarter of a million particles without compromising the computation of the gravitational 
force. Here the multifractal properties of a group of these models that incorporate different features 
of the dynamical equations governing the evolution of a matter dominated universe are compared. 
The results share important similarities with galaxy observations, such as hierarchical clustering and 
apparent bifractal geometry. They also provide insights concerning possible constraints on length 
and time scales for fractal structure. They clearly demonstrate that fractal geometry evolves in the 
fj, (position, velocity) space. The observed properties are simply a shadow (projection) of higher 
dimensional structure. 
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I. INTRODUCTION 
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Starting with the work of Vaucoulours,[l| questions 
have been posed about the geometric properties of the 
distribution of matter in the universe. Necessary for 
all modern cosmologies are the assumptions of homo- 
geneity and isotropy of the mass distribution on large 
length scales. [2] However, observations have shown the 
existence of very large structures such as super-clusters 
and voids. [3] Moreover, as technology has advanced, so 
has the length scale of the largest observed structures. 
(For a recent review see Jones et. al.Q.) It was proposed 
by Mandelbrot, based on work of Peebles, 0] that the 
matter distribution in the universe is fractal. Support for 
this conjecture came primarily from the computation of 
the pair correlation function for the positions of galaxies, 
as well as from direct observation. The fact that the cor- 
relation function was well represented by a power law in 
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the intergalactic distance seemed to support the fractal 
conjecture. In fact, in the past it has been argued by 
Pietronero and colleagues that the universe may even be 
fractal on all scales. Of course, if this were true it 
would wreak havoc with the conclusions of cosmological 
theory. 

McCauley has looked at this issue from a few different 
perspectives. He points out that since the power law be- 
havior of the correlation function is only quantitatively 
correct over a finite length scale, the universe could not 
be a simple (mono) fractal. Q Then the logical next ques- 
tion is whether or not the matter distribution is multi- 
fractal. By considering counts in cells, Bailin and Scha- 
effer hypothesized some time ago that the distribution of 
matter is approximately bifractal, i.e. a superposition of 
fractals with two distinct scaling laws.[j| McCauley has 
also addressed the possibility of inhomogeneity. Work- 
ing with Martin Kerscher, he investigated whether the 
universe has multifractal geometry. His approach was to 
examine the point-wise dimension [i| by looking for lo- 
cal scaling of the density around individual galaxies in 
two catalogues. Their conclusion was that, even if local 
scaling were present, there are large fluctuations in the 
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scaling law. Moreover, the sample sizes were not suffi- 
cient to be able to extract good local scaling exponents. 
Although larger galaxy catalogues have become available 
since their analysis, they do not meet the restrictive cri- 
teria of sufficient size, as well as uniform extension in 
all directions, necessary to measure local dimensions. [3] 
Their final conclusion was that, while the geometry of the 
observed universe is certainly not monofractal, it was not 
possible to irrefutably conclude whether it is multifractal, 
or whether the assumptions of homogeneity and isotropy 
prevail at large length scales. 

If current observations aren't able to let us determine 
the geometry of the universe, then we need to turn to 
theory and simulation. For the most part, we are con- 
cerned with the evolution of matter after the period of 
recombination. Consequently the Hubble expansion has 
slowed down sufficiently that Newtonian dynamics is ap- 
plicable, at least within a finite region of space. A 
number of theoretical approaches to the computation of 
fractal dimensions have been investigated. Each of them 
is predicated on some external assumption which is not 
yet verified. For example, de Vega et. al. assume that 
the universe is close to a thermodynamic critical point 
[ll] | , while Grujic explores a field theory where vacuum 
energy predominates over inert matter, and the latter is 
assumed to have a fractal structuring [13]. An exam- 
ination of the recent literature reveals that theory has 
not converged on a compelling and uniformly accepted 
theory of fractal structure in the universe. 

In the last few decades dynamical N-body simula- 
tion of cold dark matter CDM has experienced rapid 
advances due to improvements in both algorithms and 
technology. [3, [3] It is now possible to carry out grav- 
itational N-body simulations with upwards of 10 9 point 
mass particles. However, in order to employ simulation 
methods for systems evolving over cosmological time, 
it is necessary to compromise the representation of the 
gravitational interaction over both long and short length 
scales. For example, tree methods are frequently em- 
ployed for large separations, typically the gravitational 
field is computed from a grid, and a short range cut-off is 
employed to control the singularity in the Newtonian pair 
potential. [3, [3 Unfortunately, even if the simulations 
were perfect, a system of even 10 9 particles provides only 
10 3 particles per dimension and would thus be insufficient 
to investigate the fractal geometry with confidence. 

As a logical consequence of these difficulties it was 
natural that physicists would look to lower dimensional 
models for insight. Although this sacrifices the correct 
dynamics, it provides an arena where accurate compu- 
tations with large numbers of particles can be carried 
out for significant cosmological time. It is hoped that 
insights gained from making this trade-off are beneficial. 
In one dimension, Newtonian gravity corresponds to a 
system of infinitesimally thin, parallel, mass sheets of in- 
finite spatial extent. Since there is no curvature in a 
1+1 dimensional gravitational system, we cannot expect 
to obtain equations of motion directly from general rela- 



tivity. [15] Then a question arises concerning the inclu- 
sion of the Hubble flow into the dynamical formulation. 
This has been addressed in two ways: In the earliest, 
carried out by Rouet and Feix, the scale function was di- 
rectly inserted into the one dimensional dynamics. [iH fl7| 
Alternatively, starting with the usual three dimensional 
equations of motion and embedding a stratified mass dis- 
tribution, Fanelli and Aurell obtained a similar set of 
coupled differential equations for the evolution of the sys- 
tem in phase space. [la] While the approaches are differ- 
ent, from the standpoint of mathematics the two models 
are very similar and differ only in the values of a single 
fixed parameter, the effective friction constant. Follow- 
ing Fanelli, [T§] we refer to the former as the RF model 
and the latter as the Quintic, or Q, model. 

By introducing scaling in both position and time, in 
each model autonomous equations of motion are obtained 
in the comoving frame. In addition to the contribution 
for the gravitational field, there is now a background 
term, corresponding to a constant negative mass distri- 
bution, and a linear friction term. By eliminating the 
friction term a Hamiltonian version can also be con- 
structed. At least three other one dimensional mod- 
els have also been investigated, one consisting of New- 
tonian mass sheets which stick together whenever they 
cross, [3] one evolved by directly integrating the Zel- 
dovich equations, [13, HH and the continuous system sat- 
isfying Burger's equation [13] ■ In addition, fractal behav- 
ior has been studied in the autonomous one dimensional 
system where there is no background Hubble flow.j23l.l24j 

In dynamical simulations, both the RF and Quintic 
model clearly manifest the development of hierarchical 
clustering in both configuration and \i space (the projec- 
tion of the phase space on the position- velocity plane) . In 
common with the observation of galaxy positions, as time 
evolves both dense clusters and relatively empty regions 
(voids) develop. In their seminal work, by computing the 
box counting dimension for the RF model, Rouet and 
Feix were able to directly demonstrate the formation of 
fractal structure. {TEES] They found a value of about 0.6 
for the box counting dimension of the well evolved mass 
points in the configuration space, indicating the forma- 
tion of a robust fractal geometry. In a later work, Miller 
and Rouet investigated the generalized dimension of the 
RF model. [13] In common with the analysis of galaxy ob- 
servations by Bailin and Schaeffer Q they found evidence 
for bi-fractal geometry. Although they did not compute 
actual dimensions, later Fanelli et. al. also found a sug- 
gestion of bi-fractal behavior in the model without fric- 
tion [26j. It is then not surprising that the autonomous, 
isolated, gravitational system which does not incorporate 
the Hubble flow also manifests fractal behavior for short 
times as long as the virial ratio realized by the initial 
conditions is very small. [23. Hij In addition, in the one 
dimensional model of turbulence governed by Burger's 
equation, the formation of shocks has the appearance of 
density singularities that are similar to the clusters found 
in the RF and Quintic model. A type of bi-fractal geom- 
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etry has also been demonstrated for this system. (22| 

Below we present the results of our recent investigation 
of multifractal properties of the Quintic model and the 
model without friction. In section II we will first describe 
the systems. Then we will give a straightforward deriva- 
tion of the equations of motion and explain how they 
differ from the other models mentioned above. In section 
III we will explain how the simulations were carried out 
and describe their qualitative features. In section IV we 
define the generalized dimension and other fractal mea- 
sures and present our approach for computing them. In 
section V we will present the results of the multifractal 
analysis. Finally conclusions will be presented in section 
VI. 



II. DESCRIPTION OF THE SYSTEM 

We consider a one dimensional gravitational system 
(OGS) of TV parallel, infinitesimally thin, mass sheets 
each with mass per unit area m . The sheets are con- 
strained to move perpendicular to their surface. There- 
fore we can construct a Cartesian axis x perpendicular 
to the N sheets. We imagine that the sheets are labeled 
by j = 1,...,N. From its intersection with the x-axis, 
each sheet can be assigned the coordinate Xj and the 
corresponding velocity Vj. From Gauss' law, the force on 
each sheet is a constant proportional to the net difference 
between the total "mass" (really surface density) to its 
right and left. The equations of motion are 

^■=Vi, i£ = E i = 2nmG[N R , i -N L , i \ (1) 

where Ei is the gravitational field experienced by the 
i th particle and N^i (JVi,») is the number of particles 
(sheets) to its right (left) . Because the acceleration of 
each particle is constant between crossings, the equations 
of motion can be integrated exactly. Therefore it is not 
necessary to use numerical integration to follow the tra- 
jectory in phase space. As a consequence it was possible 
to study the system dynamics on the earliest computers 
and it can be considered the gravitational analogue of 
the Fermi-Pasta-Ulam system. [27j It was first employed 
by Lecar and Cohen to investigate the relaxation of an 
N body gravitational system. [28| Although it was first 
thought that the TV body system would reach equilib- 
rium with a relaxation time proportional to N 2 ,[29( this 
was not born out by simulations. [3p| Partially because 
of its reluctance to reach equilibrium, both single and 
two component versions of the system have been studied 
exhaustively in recent years. [H [13 (H [H] HI] [H] Most 
recently it has been demonstrated that, for short times 
and special initial conditions, the system evolution can 
be modeled by an exactly integrable system. [33] 

To construct and explore a cosmological version of the 
OGS we introduce the scale factor A(t) for a matter dom- 
inated universe. Q Moreover, as discussed in the intro- 



duction, we are interested in the development of den- 
sity fluctuations following the time of recombination so 
that electromagnetic forces can be ignored. For that, 
and later, epochs, the Hubble expansion has slowed suf- 
ficiently that Newtonian dynamics provides an adequate 
representation of the motion in a finite region. [l(| Then, 
in a 3 + 1 dimensional universe, the Newtonian equations 
governing a mass point are simply 

dr dv „ . . 

where, here, E 9 (r,t) is the gravitational field. We wish to 
follow the motion in a frame of reference where the aver- 
age density remains constant, i.e. the co-moving frame. 
Therefore we transform to a new space coordinate which 
scales the distance according to A(t). Writing r = A(t)x 
we obtain 

d 2 x 2 dA dx 1 d 2 A 1 . . 

^ + AMM + A^~^ t) (3) 

where, in the above we have taken advantage of the in- 
verse square dependence of the gravitational field to write 
E g (x, t) = ^p-E 9 (r,i) where the functional dependence 
is preserved. In a matter dominated (Einstein deSitter) 
universe we find that 



A(t) = (J-y , Pb (t) = (GnGt 2 )- 1 (4) 

where to is some arbitrary initial time corresponding, say, 
to the epoch of recombination, G is the universal gravita- 
tional constant, and pb(t) is the average, uniform, density 
frequently referred to as the background density. These 
results can be obtained directly from Eq|3]by noting that 
if the density is uniform, so that all matter is moving 
with the Hubble flow, the first two terms in Eq|3] van- 
ish whereas the third term (times A) must be equated to 
the gravitational field resulting from the uniformly dis- 
tributed mass contained within a sphere of radius A\x.\. 
Then the third term of Eq|3] is simply the contribution 
arising by subtracting the field due to the background 
density from the sphere. [1] Noting that A 3 pb(t) = pb(to) 
forces the result. Alternatively, also for the case of uni- 
form density, taking the divergence of each side of Eq. [3] 
and asserting the Poisson equation forces the same result. 
Thus the Friedman scaling is consistent with the coupling 
of a uniform Hubble flow with Newtonian dynamics. [3| 

For computational purposes it is useful to obtain au- 
tonomous equations of motion which do not depend ex- 
plicitly on the time. This can be effectively accomplished 
[La . LL3] by transforming the time coordinate according to 

dt = B(t)dr, B(t) = — (5) 

to 
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yielding the autonomous equations 

d 2 x 1 rfx 2 „ , , 

^ + ^-9T = E * (x) - (6) 

For the Newtonian dynamics considered here it is nec- 
essary to confine our attention to a bounded region of 
space, ft. It is customary to choose a cube for fl and as- 
sume periodic boundary conditions . Thus our equations 
correspond to a dissipative dynamical system in the co- 
moving frame with friction constant 1/Mq and with forces 
arising from fluctuations in the local density with respect 
to a uniform background. 

For the special case of a stratified mass distribution 
the local density at time t is given by 

p{x,t)=J2m j (t)6(x-x j ) (7) 

where rrij(t) is the mass per unit area of the j th sheet 
and, from symmetry, the gravitational field only has a 
component in the x direction. Then, for the special case 
of equal masses irij(t) = m(t), the correct form of the 
gravitational field occurring on the right hand side of Eq 
[6] at the location of particle i is then 



where a 2 is the variance of the velocity at the initial time, 
Vt = &v = a/V3 is the thermal velocity, and a is the 
maximal absolute velocity value given to a particle when 
the initial distribution of velocities is uniform on [—a, a]. 
Then, with the further requirement that L — nXj, < 
?i < N, in these units our equations are 

There is still a final issue that we have to address. If we 
assume for the moment that the sheets are uniformly dis- 
tributed and moving with the Hubble flow, then the first 
two terms above are zero. Unfortunately, if we take the 
average over each remaining term, we find that they dif- 
fer by a factor of three. This seeming discrepancy arises 
because we started with spherical symmetry about an ar- 
bitrary point for the Hubble flow and are now imposing 
axial symmetry. If we imagine that we are in a local re- 
gion with a stratified geometry, the symmetry is different 
and this has to be reflected in the background term. This 
apparent discrepancy can be rectified by multiplying the 
term in Xi by three. Our final equations of evolution are 
then 



E g (xi) = 2irm(t )G[N R , l - N L>i ] 



(8) 



since we already implicitly accounted for the fact that 
rrij{t) = m(t[))/A 2 in EqEl The equations are further 
simplified by establishing the connection between m(rj ) 
and the background density at the initial time, pb(to). Let 
us assume that we have 2N particles (sheets) confined 
within a slab with width 2L, i.e. —2L < x < 2L. Then 

p b (t ) = (GnGtiy 1 = (j-^j rofo), (9) 
we may express the field by 

2 (L 



Eg{ Xi ) = 2-Km{t )G[N R ,i-NL,i] = ^ ) [Nn,i-N L ,i] 

(10) 

and the equations of motion for a particle in the system 
now read 



dr 2 Mo dr 9t 2 3t% \N 



[NR,i-N L ,i]. (11) 



It is convenient to refer to Jeans theory for the final choice 
of units of time and length [38] 



— <r„t n , (12) 



d Xi 
~d^ 



1 dxi 
V^~d^ 



iV 



L,i\ 



(14) 



The description is completed by assuming that the sys- 
tem satisfies periodic boundary conditions on the interval 
2L, i.e. when a particle leaves the primitive cell defined 
by — 2L < x < 2L on the right, it re-enters at the left 
hand boundary with the identical velocity. Note that, in 
computing the field, we do not attempt to include con- 
tributions from the images of Xi outside of the primitive 
cell. 

Finally, we mention that the RF model is obtained 
from the reverse sequence where one first restricts the ge- 
ometry to 1+1 dimensions and then introduces the trans- 
formation to the comoving frame. In this approach it is 
not necessary to make the adjustment in the coefficient 
of Xi as we did here. This is quickly seen by noting that 
the divergence of x is three times greater than the diver- 
gence of xx which one would obtain by directly starting 
with the one dimensional model. However, in the RF 
model, the coefficient of the first derivative term (the 
friction constant) is 1/ V2 instead of 1/ y/E. This simply 
illustrates that, since there is no curvature in a 1+1 di- 
mensional universe, there is a degree of arbitrariness in 
choosing the final model. It cannot be obtained solely 
from General Relati vity . For a discussion of this point 
see Mann et. al. . [15J A Hamiltonian version can also 
be considered by setting the friction constant to zero. 
In their earlier work, using the linearized Vlasov-Poisson 
equations, Rouet and Feix carried out a stability analy- 
sis of the model without friction. They determined that 
the system followed the expected behavior, i.e. when the 
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system size is greater than the Jean's length, instability 
occurs and clustering becomes possible. LL3] We men- 
tion that, when the friction term is not present, both the 
Q and RF models are identical so, with the assumption 
that the friction term won't have a large influence on 
short time, linear, stability, the analysis of the Hamilto- 
nian version applies equally to each version. 

The Vlasov-Poisson limit for the system is obtained 
by letting TV — > oo and m — > while constraining the 
density and energy (at a given time). Then, in the co- 
moving frame, the system is represented by a fluid in the 
/i space. Let f(x,v;t) represent the normalized distribu- 
tion of mass in the (x, v) phase plane at time t. From 
mass conservation, / satisfies the continuity equation 



df 

at 



df , daf 



da 



dv 



0. 



-7U 



d(j>r{x) 
dx 



(15) 



which can also be derived using the identical scaling as 
aboveS. In Eq. ([IB]). a(x,v) is the local acceleration and 
4>t is the potential function induced by the total density, 
including the effective negative background density, —pb- 
Note that both a(x,v) and <f>T are linear functional of 
f(x,v;t). Coupling Eq. (J5|) with the Poisson equation 
yields the complete Vlasov-Poisson system governing the 
evolution of f(x,v;t). Depending on the final choice of 
the friction constant, 7, the continuum limit of either the 
Quintic or RF model can be represented by EqUH Note 
that an exact solution is 



f(x,v,t) = 



Pb\ 



V2? 



■. exp 



~2cj2' = CT o cxp -7i 



(16) 

Using dynamical simulation, we will see that it is ex- 
tremely unstable when the system size exceeds the Jeans' 
length at the initial time. 

Useful information can be had without constructing 
an explicit solution. For example, we quickly find that 
the system energy decreases at a rate proportional to the 
kinetic energy. The Tsallis entropy is defined by 



■So 



1 - / f q dxdv 



fdxdv = 1 



(17) 



In the limit q — > 1, S q reduces to the usual Gibbs entropy, 
5*1 = — J J f 'In fdxdv. By asserting the Vlasov-Poisson 
evolution, we find that the Bolztmann-Gibbs entropy, Si, 
decreases at the constant rate —2j while, for q > 1, S q 
decreases exponentially in time. This tells us that the 
mass is being concentrated in regions of decreasing area 
of the phase plane, suggesting the development of struc- 
ture. These properties are immediately evident for the 
trivial solution given above. By imposing a Euclidean 
metric in the phase plane, we can also investigate local 
properties such as the directions of maximum stretching 
and compression, as well as the local vorticity. We eas- 
ily find that the rate of separation between two nearby 
points is a maximum in the direction given by 



tan(20) = (l + p + p h )h, 



(18) 



where 9 defines the angle with the coordinate axis in the 
p space. Thus, in regions of low density, we expect to see 
lines of mass being stretched with constant positive slope 
in the phase plane. We will see below that this prediction 
is accurately born out by simulations. 



III. SIMULATIONS 

An attraction of these one-dimensional gravitational 
systems is their ease of simulation. In both the au- 
tonomous (purely Newtonian without expansion) and RF 
models it is possible to integrate the motion of the in- 
dividual particles between crossings analytically. Then 
the temporal evolution of the system can be obtained by 
following the successive crossings of the individual, adja- 
cent, particle trajectories. This is true as well for the Q 
model. If we let tji = Xj+i — x,, where we have ordered 
the particle labels from left to right, then we find that 
the differential equation for each yi is the same, namely 



dT 2 



1 dy,j 
V6 dr 



2n 
' TV' 



(19) 



The general solution of the homogenious version of 
E. [19] is a sum of exponentials. By including the par- 
ticular solution of the inhomogenius equation (simply a 
constant) we obtain a fifth order algebraic equation in 
u = cxp(r/\/6)- Hence the name Q, or Quintic, model. 
These can be solved numerically in terms of the initial 
conditions by analytically bounding the roots and em- 
ploying the Newton- Raphson method. (3^1 (Note that for 
the RF model a cubic equation is obtained.) A sophis- 
ticated, event driven, algorithm was designed to execute 
the simulations. Two important features of the algorithm 
are that it only updates the positions and velocities of a 
pair of particles when they actually cross, and it main- 
tains the correct ordering of each particle's position on 
the line. Using this algorithm we were able to carry out 
runs for significant cosmological time with large numbers 
of particles. 

Typical numerical simulations were carried out for sys- 
tems with up to TV = 2 18 particles. Initial conditions were 
chosen by equally spacing the particles on the line and 
randomly choosing their velocities from a uniform dis- 
tribution within a fixed interval. For historical reasons 
we call this a waterbag. Other initial conditions, such 
as Normally distributed velocities, as well as a Brown- 
ian motion representation, which more closely imitates 
a cosmological setting, were also investigated. The sim- 
ulations were carried forward for approximately fifteen 
dimensionless time units. 

In Fig. [TJ we present a visualization of a typical run 
with 2 17 particles. The system consists of the Q model 



with an initial waterbag distribution in the p space. Ini- 
tially the velocity spread is (-12.5, 12.5) in the dimen- 
sionless units employed here and the system contains ap- 
proximately 16,000 Jeans' lengths. In the left column 
we present a histogram of the particles positions at in- 
creasing time frames, while on the right we display the 
corresponding particle locations in \i (position, velocity) 
space. It is clear from the panels that hierarchical clus- 
tering is occurring, i.e. small clusters are joining together 
to form larger ones, so the clustering mechanism is "bot- 
toms up" [2]. The first clusters are roughly the size of a 
Jean's length and seem to appear at about r = 6 and 
there are many while by t = 14 there are on the order 
of 30 clusters. In the /x space we observe that between 
the clusters matter is distributed along linear paths. As 
time progresses the size of the linear segments increases. 
The behavior of these under-dense regions is governed by 
the stretching in fj, space predicted by Vlasov theory ex- 
plained above (see before Eq.JTH}). The slopes of the seg- 
ments are the same and in quantitative agreement with 
Eq. (fl8|) . Qualitatively similar histories are obtained 
for the RF model and the model without friction (see 
below), as well as for the different boundary conditions 
mentioned above. However, there are some subtle dif- 
ferences. In Fig. [2] we zoom in and show a sequence of 
magnified inserts from the ^space at time T=14. The 
hierarchical structure observed in these models suggests 
the existence of a fractal structure, but careful analysis is 
required to determine if this is correct. Simulations have 
also been performed where the system size is less than the 
Jean's length. The results support the standard stability 
analysis in that hierarchical clustering is not observed for 
these initial conditions. 

Historically, power law behavior in the density-density 
correlation funtion has been taken as the most important 
signature of self similar behavior of the distribution of 
galaxy positions. [3] In Figure [3] we provide a log-log plot 
of the correlation function C(r) at T=14 defined by 



the multifractal analysis described table U below in some 
detail. 



IV. FRACTAL MEASURES 

It is natural to assume that the apparently self similar 
structure which develops in the phase plane (see Figs Jl|2|) 
as time evolves has fractal geometry, but we will see that 
things aren't so simple. In their earlier study of the RF 
model, Rouet and Feix found a box counting dimension 
for the particle positions of about 0.6 for an initial wa- 
terbag distribution (uniform on a rectangle in the phase 
plane-see above) and a fractal dimension of about 0.8 in 
the configuration space (i.e. of the projection of the set of 
points in \x space on the position axis). fl7l |. As far as we 
know, Bailin and Shaeffer were the first to suggest that 
the distribution of galaxy positions are consistent with a 
bifractal geometry [40j . Their idea was that the geome- 
try of the galaxy distribution was different in the clusters 
and voids and, as a first approximation, this could be rep- 
resented as a superposition of two independent fractals. 
Of course, their analysis was restricted solely to galaxy 
positions. Since the structures which evolve are strongly 
inhomogeneous, to gain further insight we have carried 
out a multi-fractal analysis [Il| in both the phase plane 
and the position coordinate. 

The multifractal formalism shares a number of features 
with thermodynamics. To implement it we partitioned 
each space configuration space and /i space, into cells of 
length I. At each time of observation in the simulation, a 
measure m = Ni(t)/N was assigned to cell i, where Ni(t) 
is the population of cell i at time t and N is the total 
number of particles in the simulation. The generalized 
dimension of order q is defined by [4l| 



D„ 



-lim^. 
1 i^o lal 



C q = Xtf. 



(20) 



C{r) =< Sp(y+r)6p(y) >~ — 



r+A 



where particles i and j are such that L/4 < Xi,j < 3L/4 
to avoid boundary effects and A is the bin size. Note the 
existence of a scaling region from about 0.1 < ln(Z) < 
30, a range of about 2.5 decades in I. Note also the 
noise present at larger scales. Since the fluctuations are 
vanishingly small at scales on the order of the system 
size, this is most likely attributable to the presence of 
shot noise. This could be reduced by taking an ensemble 
average over many runs. By computing the slope, say 7 
, of the log-log plot in Fig. [3] we are able to obtain an 
estimate of the correlation dimension D2 = 1— 7 for a one 
dimensional system. We find that 7 = .42 for a scaling 
region of about two decades in I . This suggests that 
the correlation dimension is approximately 0.6, which is 
in agreement within the standard numerical error with 



?)jjs the effective partition function Do is 
yj ) )the 1 box" count : ngrdimension, D\, obtained by taking the 
limit q-L— ► 1, is the information dimension, and Z?2 is 
the correlation dimension [Il]][:9]. As q increases above 0, 
the D q provide information on the geometry of cells with 
higher population. 

In practice, it is not possible to take the limit I — > 
with a finite sample. Instead, one looks for a scaling rela- 
tion over a substantial range of In I with the expectation 
that a linear relation between lnC q and \al occurs, sug- 
gesting power law dependence of C q on I. Then, in the 
most favorable case, the slope of the linear region should 
provide the correct power and, after dividing by q — 1, 
the generalized dimension D q . As a rule, or guide, if scal- 
ing can be found either from observation or computation 
over three decades of /, then we typically infer that there 
is good evidence of fractal structure. pj Also of interest 
is the global scaling index r q , where C q ^l Tq for small I . 
It can be shown that r q and D q are related to each other 
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Figure 1: Evolution in configuration and /x space for the quintic model with 2 17 particles from T=0 to T=14. The initial 
distribution is a waterbag with velocities in the range (-12.5, 12.5) and a size of about 16,000 Jeans' lengths. 



T q for 



through a Legendre transformation by D q (l — q) 
q 7^ 1 [9(| . Here we present the results of our fractal anal 
ysis of the particle positions on the line (position only) 
and the plane (position, velocity or \i space). 

If it exists, a scaling range of I is defined as the in- 
terval on which plots of In C q versus In I are linear. Of 
course, for the special case of q = 1 we plot — E^jZn(/ii) 
vs. ln(l) to obtain the information dimension. [3| If a 
scaling range can be found, D q is obtained by taking the 
appropriate derivative. It is well established by proof and 
example that, for a normal, homogeneous, fractal, all of 
the generalized dimensions are equal, while for an inho- 
mogeneous fractal, e.g. the Henon attractor, D q+ i < D q 
[4l| . In the limit of small /, the partition function, C q (l), 
can also be decomposed into a sum of contributions from 
regions of the inhomogeneous fractal sharing the same 
point-wise dimension, a, 



C q {l) = / danp(a) exp [-/ (a)] 



where 

(HE 



(21) 



[a) is the fractal dimension of its support 
. Then if, for a range of q, a single region is 



dominant, we find a simple relation between the global 
index, r g , and a, 



T q = aq- f(a). 



(22) 



Recall that a Euclidean metric is imposed on the /it- 
space. Because, in our units, uij = 1, we divide it 
into cells such that Ax = Av. Then, as L is large 
(L fa 10,000), so is the extent of the partition in the 
velocity space. On this scale the initial distribution of 
particules appears as a line so that the initial dimen- 
sion is also about unity. In fact, initially the velocity of 
the particles is a perturbation. It is not large enough 
to allow particles to cross the entire system in a unit of 
time. During the initial period the virial ratio rapidly 
oscillates. After a while the granularity of the system de- 
stroys the approximate symmetry of the initial /z-space 
distribution. Breaking the symmetry leads to the short 
time dissipative mixing that results in the separation of 
the system into clusters. Although the embedding dimen- 
sion is different, the behavior of the distribution of points 
in configuration space is similar. The initial dimension is 



8 



5000 



> 2000 



I- 














"i . 













-1000 

56000 57000 58000 59000 60000 61000 62000 63000 64000 

(a) x 



2500 



1000 - 



58450 58500 58550 58600 58650 58700 58750 58800 58850 

(b) X 

Figure 2: Consecutive expansions (zooms) on successive small 
squares selected in the /i space panels at T=14 for the quintic 
model. They have the appearance of a random fractal which 
suggests self-similarity. 
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Figure 3: Figure 3. The correlation function at T=14 which 
shows a scaling region from about 0.3 to about 30, a range of 
about 2 decades with a scaling exponent 7 = .42. 



nearly one until clustering commences. At this time the 
dimensions in /i-space and configuration space separate. 



As time progresses, however, for the initial conditions 
discussed above, for q > typically two independent scal- 
ing regimes developed. Of course this is in addition to 
the trivial scaling regions obtained for very small / , cor- 
responding to isolated points, and to large I on the or- 
der of the system size, for which the matter distribution 
looks smooth. The observed size of each scaling range de- 
pended on both the elapsed time into the simulation and 
the value of q. While the length of each scaling regime 
varied with both q and time, in some instances it was 
possible to find good scaling over up to four decades in l\ 

In Fig. 2] we provide plots of -^jlnC q versus Inl in 
/i-space for four different values of q covering most of 
the range we investigated (—5 < q < 10). To guarantee 
that the fractal structure was fully developed, we chose 
T = 14 for the time of observation and the initial con- 
ditions are those given above. For q = — 5 and q = 
we clearly observe a single, large, scaling range where 
0.5 < Inl < 8 , corresponding to about three decades in 
/. In the remaining panels (c, d), where we increase q 
to 5 , and then 10, we see a dramatic change. The large 
scaling range has split into two smaller regions separated 
at about I = 1.5, and the slope of the region with larger 
I has decreased compared with the scaling region on its 
left (—1.5 < In / < 1.5). The scaling range with larger 
1.5 < \nl < 8.5, corresponds to just over three decades 
in I. Note that in panels c and d of Fig. IH the scaling 
range with larger I is more robust. 

Now that we have identified the important scales, we 
are able to compute the generalized dimension, D q , and 
the global index, r q . In Fig. 5(a) we plot D q vs q cal- 
culated in the /i-space for the quintic model at T = 14 
for q > for the case of the larger scaling range. As 
expected, it is a decreasing function and clearly demon- 
strates multifractal behavior. Although \i space is two 
dimensional, Do is about 0.9. While Dq > D\ is typ- 
ical for multifractals (information dimension covers the 
important region), the fact that D q is strictly decreasing 
suggests greater "fractality" from increasingly overdense 
regions, i.e. the system is strongly inhomogeneous. In 
Fig. 5b we plot r q versus q in mu space. Two linear re- 
gions (0 < q < 1 and 1 < q < 10) can be distinguished 
suggesting bifractal geometry, i.e. a superposition of two 
fractals with unique values of a and f(a). 

We have seen that the one dimensional gravitational 
system reveals fractal geometry in the higher dimensional 
/i space (frequently referred to as phase space in the as- 
tronomical literature). Historically self similar behavior 
was first inferred in gravitational systems by studying 
the distribution of galaxy locations, i.e in the system's 
configuration space, and searching for powerlaw behav- 
ior in the correlation function (see above). [53] In common 
with this approach, we project the mu space distribution 
on the (position) line to study the geometry of the con- 
figuration space. In Fig. [6] a,b,c,d, we provide plots of 
l/(g— l)ln(I q ) vs ln(l) in x (configuration) space for the 
quintic model prepared as above at T=14. In common 
with the /i-space distribution, for q = —5 and q — 0, there 
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Figure 4: Scaling behavior in ^t-space. Plots of ^-j-lnC 9 ver- 
sus InZ are provided for four values of q: a) q=-5, b) q=0. c) 
q=5, d) q= 10. 



is a single, large scaling region (—3 < ln(l) < 3 for q = — 5 
and —2 < ln(l) < 6 for q = 0). Similarly, at q — 5 (c) and 
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Figure 5: Generalized dimension D q vs. q (panel a) and global 
scaling index r^vs. q (panel b) in /i-space for the quintic model 
at T=14 with q > 0. 



q = 10 (d) the existence of two scaling ranges becomes 
apparent, —4 < ln(l) < and < ln(l) < 8 for q = 5 
(about 3.3 decades), -3 < ln(l) < 1 and 1 < ln(l) < 8 for 
q = 10 (about 3 decades). Note that the scaling ranges 
are similar for both the \i and configuration space, but 
they are not identical. 

In Figure[7|we illustrate the behavior of the generalized 
dimension D q &nd global index r q for the configuration 
space of the quintic model at T = 14. Although now the 
embedding dimension is d = 1, D is about 0.7 so the 
distribution is definitely fractal. As anticipated from the 
/i-space distribution, D q is a decreasing function of q so 
it is also inhomogeneous and multifractal. In Fig. [7t> we 
plot T q vs q in configuration space, for the same system. 
In contrast with the /i-space index, here we observe a 
nearly linear function with slight convexity (decreasing 
derivative with respect to q). Perhaps this results from 
the superposition of three linear regions with distinct a 
and /(a), say with < q < 1.5; 1.5 < q < 7; 7 < q < 10, 
but this cannot be inferred from the plot without further 
analysis. 

For completeness, in Fig. [8] we examine the behavior 
of D q and r q for q < 0. In general, negative q is im- 
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Figure 6: Plots of \/{q — l)ln(C q ) vs ln{l) in configuration (a: 
) space for the quintic model prepared as above at T=14 are 
provided for four values of q: a) q = —5, b) g = 0, c) q = 5, 
d)q= 10. 

portant for revealing the geometry of low density regions 
(voids). We see that D q has a very unphysical appear- 
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Figure 7: Generalized dimension D q vs. q (panel a) and global 
scaling index r q vs. q (panel b) in configuration (x)space for 
the quintic model at T=14 with q > 0.. 

ance - it is increasing. The source of the problem can be 
determined by examining the behavior of r q . It is nearly 
constant over the range — 10 < q < with a value of 
approximately r q ~ —1, implying that the underdense 
regions are dominated by a single local dimension. 

Part of our goal is to compare how fractal geometry 
arises in a family of related models. So far we have pre- 
sented results for the quintic, or Q, model. However we 
have also carried out similar studies of the Rouet-Feix 
(RF) model, the model without friction (which can be 
obtained from either of the former by nullifying the first 
derivative contribution in Eq (fl4]) . and simply an isolated 
system without friction or background. The latter is 
a purely Newtonian model without a cosmological con- 
nection. In Table U we list the important characteris- 
tic dimensions and exponents for the the quintic model 
and the model without friction for comparison. While 
there are similarities in the fractal structure of each of 
these models, they are not the same. For example in 
the quintic model the generalized fractal dimension, D q , 
is consistently smaller in each manifold than the corre- 
sponding dimension in the model without friction. It is 
noteworth that the box counting dimension in the config- 
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Figure 8: Generalized dimension D q vs. q (panel a) and global 
scaling index r q vs. q (panel b) in configuration (x')space for 
the quintic model at T=14 with q < 0.. 



urqtion space of the quintic model is 0.69 compared with 
0.90 for the model without friction. Moreover D q is de- 
creasing more rapidly in the quintic model demonstrating 
stronger inhomogeneity. 



Table I: default 
Quintic model 



RF model 



x-space /i-space 
0.69 ±0.11 0.89 ±0.06 
0.64 ± 0.06 0.84 ± 0.07 
0.62 ±0.04 0.81 ±0.05 



x-space /i-space 
0.89 ±0.02 1.14 ±0.04 
0.88 ±0.04 1.17 ±0.04 
0.85 ±0.04 1.13 ±0.05 



The differences in dynamics between the systems can 
be seen clearly in Fig. [9] where we show the evolution 
in configuration and ^-space of the system without fric- 
tion. In common with the simulations of the quintic 
model presented above, the initial conditions are a wa- 
terbag with velocities sampled independently and con- 
fined to (—12.5, 12.5) in the dimensionless units defined 
earlier. As before, there are 2 17 particles and the simu- 
lation spans 14 units of the scaled time. The initial size 



encompasses about 10,400 Jeans lengths. Compared with 
the quintic model, note that here larger clusters form ear- 
lier and have a qualitatively different shape. Due to the 
absence of "friction" in the dynamics, the velocity spread 
is larger and there are fewer clusters at each epoch. At 
T=12 and T=13 the linear structure of the underdense 
regions (voids) in fi space connecting the clusters is ap- 
parent. Note that they all share a common slope which is 
due to stretching in the phase plane. However, equation 
fl8l is no longer adequate to provide the line slope. In 
Fig. [10] we examine the consecutive expansions (zooms) 
of the large cluster on the right hand side of the \i space 
distribution at T = 14. Once again there is qualitative 
evidence of stochastic self-similarity j§] which requires 
analysis for confirmation. To save space, here we don't 
reproduce the plots of ^yln(C g ) versus In I in each space. 
Qualitatively they are similar to the quintic case, but the 
scaling ranges are less robust. 

The behavior of D q and r q are similar to that found for 
the quintic model. However, because the scaling ranges 
are less robust, there is more noise. For example, in Fig. 
fTTk.b. we plot D q and r q versus q in configuration space 
for the model without friction. Note the similarities with 
Fig. [6] for the quintic model. The slope of r q gradually 
decreases with increasing q. There appears to be three 
linear segments, from (0, 3), (3, 5), and (5,10) suggesting 
contributions from three possible fractal scaling regions. 



V. DISCUSSION OF RESULTS 

As mentioned earlier, it is well established that, for a 
regular multifractal, the generalized dimension, D a , is a 
decreasing function of its argument. [42j[9J] Therefore, for 
q < , it would be incorrect to interpret the simulation 
results as true generalized dimensions. There must be an 
alternative explanation for the behavior we observe. The 
picture for positive q is rather different. The two scal- 
ing regions give completely different results. Although 
one would suspect from the definition of the generalized 
dimension that the scaling region with smaller I would 
give the correct one, this is hard to accept since typical 
plots of the function D q are still increasing until about 
q = 2 (not shown) for this range. On the other hand, 
the second, larger, scaling region manifests a well be- 
haved decreasing function which appears to approach a 
constant value of D q = 0.63 for q < 10. 

It is interesting that we have observed similar behavior 
with a well characterized, textbook, fractal that is dis- 
cussed in numerous sources. (42| As a test of our compu- 
tational approach we simulated the multiplicative bino- 
mial process. For this multifractal r q and D q are known 
precisely. [12] When we carried out the fractal analysis 
using the methods described above, we also found two 
scaling regions. What is most striking is that the scaling 
region with larger values of I yielded a r q (and there- 
fore D q ) which agreed to within numerical error with the 
theoretical prediction! This type of behavior has been 
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Figure 9: Evolution in configuration and space for the model without friction with 2 17 particles from T=0 to T=14. The initial 
distribution is a waterbag with velocities in the range (-14, 14) in the dimensionless units employed here and a size of about 
10,400 Jeans' lengths. 



manifested in other simulations of fractal sets. The for- 
mation of the smaller scaling range has been attributed 
to the presence of noise in the data. The existence of a 
second scaling range has been rigorously shown to arise 
when Gaussian noise is added to a standard fractal. [43| 
In the simulation of the multiplicative binomial process, 
the source of the noise is the random location of the data 
points within the smallest bins. At this time the source of 
the apparent noise in the one dimensional gravitational 
simulations is not precisely known. While it may arise 
simply from numerical considerations, there are alterna- 
tive possible explanations. For example, noise may arise 
from sub-Jeans length fluctuations in the initial data, or 
from other small scale features of the initial state. In 
addition to the box-counting approach, other methods 
based on the correlation function formalism that employ 
the point-wise dimension can also be used to investigate 
the fractal properties of the system. It was also used 
in our investigation and gives similar results to the box 
counting method reported here. 



ena to the mutiplicative binomial process or the systems 
with injected noise. [43| Then how do we explain the sur- 
prising and counterintuitive results for q < 0? Since for 
negative q we obtain a nearly constant value for r q from 
each scaling region, it seems safe to assume that a region 
of the data characterized by a simple fractal behavior 
has the dominant influence. Moreover, since it only in- 
volves q < 0, it represents the regions of low density, 
i.e. the voids. Referring to the multifractal formalism, 
we pointed out earlier that in a region where a single 
structure dominates, 



aq- f(a), 



(23) 



where a is the strength of the local singularity, or the lo- 
cal pointwise dimension,and f(a) is the Hausdorf dimen- 
sion of its support. [Hill Then the computations show 
that for q < we must have a = and / = 0.9. This 
suggests that the results for negative q are dominated by 
regions of such low density that widely separated, "iso- 
lated", particles are responsible for the spurious behavior 
of D q . At this time this is simply a conjecture which 



For q > we appear to be seeing a similar phenom- needs to be investigated with further computation. 



13 




-25900 



70000 75000 60000 B500Q 

(a) s 

-2000 i , 1 1 



95000 1 j: 



-5000 






> 






-6000 


■ il?' 

\ 


f . r • 






-7000 


\ 





-8000 - 



-9000 ' 

77000 



78000 



79000 



(b) 



80000 
X 



81000 



82000 



83000 



Figure 10: Consecutive expansions (zooms) of successive 
small squares (square within a square) selected in the /i space 
panels at T=14. They also have the appearance of a random 
fractal which suggests self-similarity. 



Finally let's reconsider the plot of r q versus q (Fig. [5b) 
for the larger scaling region with q > 0. As we mentioned 
earlier, for < q < 1.5 and for 2 < q < 10 the curve 
appears linear with different slopes in each region. This 
may be the manifestation of bifractality first discussed by 
Bailin and Shaffer for galaxy positions. The first interval 
may represent the true fractal structure of the under- 
dense regions, while the dense clusters are dominant for 
the larger q values. Bifractal behavior is not unique to 
gravitational systems. It occurs in a number of well stud- 
ied model systems, for example as a superposition of two 
Cantor sets,[H] or from the truncation of Levy flights 
1451141 • 



VI. CONCLUSIONS 

We have seen that one dimensional models develop hi- 
erarchical structure and manifest robust scaling behavior 
over particular length and time scales. In addition, for 
a number of reasons, they have an enormous computa- 
tional advantage over higher dimensional models. First, 
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Figure 11: Generalized dimension D q vs. q (panel a) and 
global scaling index r 9 vs. q (panel b) in configuration (x)space 
for the quintic model at T=14 with q < 0. 



it is possible to study the evolution of large systems with 
on the order of 2 17 particles per dimension. Second, in 
contrast with three dimensional N-body simulations, the 
evolution can be followed for long times without com- 
promising the dynamics. In particular, the force is al- 
ways represented with the accuracy of the computer - 
there is no softening. As a consequence it is possible to 
investigate the prospect for fractal geometry with some 
confidence. 

In common with 3+1 dimensional cosmologies, with 
the inclusion of the Hubble expansion and the transfor- 
mation to the comoving frame, both the 1+1 dimensional 
Quintic and RF models reveal the formation of dense 
clusters and voids. They also show evidence for bifractal 
geometry. This may be a consequence of dynamical in- 
stability that results in the separation of the system into 
regions of high and low density. An interesting obser- 
vation is that the lower bound of the length scale that 
supports the trivial space dimension of unity in the con- 
figuration space grows with time. To the extent that 
similar behavior occurs in the 3+1 dimensional universe, 
this lends support for the standard cosmological model 
on sufficiently large scales. It also suggests that the scale 
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size for homogeneity will grow with time. Eventually this 
may be testable with observation. 

We have seen that the system shows evidence of two 
nontrivial scaling regions. The type of anomalous be- 
havior of D q and T q in the scaling region with a finer 
partition was also found in the standard multiplicative 
binomial process with a similar sample size. Even with 
the large sample size employed here, this would suggest 
that it is a finite size effect. In that respect it will be 
interesting to perform simulations dealing directly with 
the distribution function given by the Vlasov equation. 
Then the low density region will be well described. This 
may also be true for the region of negative q where a 
fractal analysis forces us to conclude that the point-wise 
dimension vanishes. Computations with different bound- 
ary conditions reveal similar behavior, but this work is 
only in the preliminary stages. Future work will include 
the investigation of the influence of correlation in the ini- 
tial conditions, the consideration of initial conditions of 
the type employed in current cosmological N-body simu- 
lations [H] , and the study of the fractal geometry of the 
under and over dense regions. We plan to compare out 
findings with three dimensional studies of the distribu- 
tion of voids and halos [49] . 

In addition to these important features, the one dimen- 
sional simulations unequivocally demonstrate that the es- 



sential structure formation is taking place in the /i space 
(phase space for astronomers). Thus the apparent frac- 
tal geometry that we observe in configuration space is 
simply a shadow (projection) of the higher dimensional 
structure. The development of structure in fi space and 
its apparent fractal geometry may be the most significant 
result of these studies. Since evolution in the six dimen- 
sional /i space of the observable universe is beyond our 
current capability, the study of lower dimensional models 
is an important guide for understanding the important 
features of the higher dimensional evolution. For exam- 
ple, we have shown analytically that dynamical "stretch- 
ing" of the [i space geometry is responsible for the for- 
mation of underdense regions (voids) and, consequently, 
the concentration of mass in regions of decreasing area. 
We are currently extending this analysis to the study of 
structure formation in the more realistic 3+1 dimensional 
manifold. 
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